home *** CD-ROM | disk | FTP | other *** search
/ Computer Select (Limited Edition) / Computer Select.iso / dobbs / v17n04 / wavelet.exe / WAVELET.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-17  |  7.4 KB  |  268 lines
  1. /* WAVELET.C */
  2.  
  3. #include <math.h>
  4.  
  5. typedef enum { DECOMP, RECON } wavetype;
  6.  
  7. #include "wavelet.h"
  8.  
  9. void WaveletCoeffs(double alpha, double beta, double *wavecoeffs)
  10. {
  11.     double tcosa, tcosb, tsina, tsinb;
  12.     char i;
  13.     /* precalculate cosine of alpha and sine of beta to reduce */
  14.     /* processing time */
  15.     tcosa = cos(alpha);
  16.     tcosb = cos(beta);
  17.     tsina = sin(alpha);
  18.     tsinb = sin(beta);
  19.  
  20.     /* calculate first two wavelet coefficients  a = a(-2) and b = a(-1) */
  21.     wavecoeffs[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
  22.                                                                             + 2.0 * tsinb * tcosa) / 4.0;
  23.     wavecoeffs[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
  24.                                                                             - 2.0 * tsinb * tcosa) / 4.0;
  25.  
  26.     /* precalculate cosine and sine of alpha minus beta to reduce */
  27.     /* processing time */
  28.     tcosa = cos(alpha - beta);
  29.     tsina = sin(alpha - beta);
  30.  
  31.     /* calculate last four wavelet coefficients  c = a(0), d = a(1), */
  32.     /* e = a(2), and f = a(3) */
  33.     wavecoeffs[2]  = (1.0 + tcosa + tsina) / 2.0;
  34.     wavecoeffs[3]  = (1.0 + tcosa - tsina) / 2.0;
  35.     wavecoeffs[4]  = 1 - wavecoeffs[0] - wavecoeffs[2];
  36.     wavecoeffs[5]  = 1 - wavecoeffs[1] - wavecoeffs[3];
  37.  
  38.     /* zero out very small coefficient values caused by truncation error */
  39.     for (i = 0; i < 6; i++)
  40.     {
  41.         if (fabs(wavecoeffs[i]) < 1.0e-15)
  42.             wavecoeffs[i] = 0.0;
  43.     }
  44.  
  45. }
  46.  
  47.  
  48. char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
  49.                                 double *Hfilter, wavetype transform)
  50. {
  51.     char i, j, k, filterlength;
  52.  
  53.     /* find the first non-zero wavelet coefficient */
  54.     i = 0;
  55.     while(wavecoeffs[i] == 0.0)
  56.         i++;
  57.  
  58.     /* find the last non-zero wavelet coefficient */
  59.     j = 5;
  60.     while(wavecoeffs[j] == 0.0)
  61.         j--;
  62.  
  63.     /* form the decomposition filters h~ and g~ or the reconstruction
  64.          filters h and g.  the division by 2 in the construction
  65.          of the decomposition filters is for normalization */
  66.     filterlength = j - i + 1;
  67.     for(k = 0; k < filterlength; k++)
  68.     {
  69.         if (transform == DECOMP)
  70.         {
  71.             Lfilter[k] = wavecoeffs[j--] / 2.0;
  72.             Hfilter[k] = (double) (((i++ & 0x01) * 2) - 1) * wavecoeffs[i] / 2.0;
  73.         }
  74.         else
  75.         {
  76.             Lfilter[k] = wavecoeffs[i++];
  77.             Hfilter[k] = (double) (((j-- & 0x01) * 2) - 1) * wavecoeffs[j];
  78.         }
  79.     }
  80.  
  81.     /* clear out the additional array locations, if any */
  82.     while (k < 6)
  83.     {
  84.         Lfilter[k] = 0.0;
  85.         Hfilter[k++] = 0.0;
  86.     }
  87.  
  88.     return filterlength;
  89. }
  90.  
  91.  
  92. double DotP(double *data, double *filter, char filtlen)
  93. {
  94.     char i;
  95.     double sum;
  96.  
  97.     sum = 0.0;
  98.     for (i = 0; i < filtlen; i++)
  99.         sum += *data-- * *filter++; /* decreasing data pointer is moving */
  100.                                                                 /* backward in time */
  101.     return sum;
  102. }
  103.  
  104.  
  105. void ConvolveDec2(double *input_sequence, int inp_length,
  106.                     double *filter, char filtlen, double *output_sequence)
  107.     /* convolve the input sequence with the filter and decimate by two */
  108. {
  109.     int i;
  110.     char shortlen, offset;
  111.  
  112.     for(i = 0; (i <= inp_length + 8) && ((i - filtlen) <= inp_length + 8); i += 2)
  113.     {
  114.         if (i < filtlen)
  115.             *output_sequence++ = DotP(input_sequence + i, filter, i + 1);
  116.         else if (i > (inp_length + 5))
  117.         {
  118.             shortlen = filtlen - (char) (i - inp_length - 4);
  119.             offset = (char) (i - inp_length - 4);
  120.             *output_sequence++ = DotP(input_sequence + inp_length + 4,
  121.                                                                                 filter + offset, shortlen);
  122.         }
  123.         else
  124.             *output_sequence++ = DotP(input_sequence + i, filter, filtlen);
  125.     }
  126.  
  127. }
  128.  
  129.  
  130. int DecomposeBranches(double *In, int Inlen, double *Lfilter,
  131.         double *Hfilter, char filtlen, double *OutL, double *OutH)
  132.     /* take input data and filters and form two branches of have the
  133.        original length. length of branches is returned */
  134. {
  135.     ConvolveDec2(In, Inlen, Lfilter, filtlen, OutL);
  136.     ConvolveDec2(In, Inlen, Hfilter, filtlen, OutH);
  137.     return (Inlen / 2);
  138. }
  139.  
  140.  
  141. void WaveletDecomposition(double *InData, int Inlength, double *Lfilter,
  142.             double *Hfilter, char filtlen, char levels, double **OutData)
  143.     /* assumes that the input data has 2 ^ (levels) data points per unit
  144.          interval.  first InData is input signal; all others are the
  145.          intermediate approximation coefficients */
  146. {
  147.     char i;
  148.  
  149.     for (i = 0; i < levels; i++)
  150.     {
  151.         Inlength = DecomposeBranches(InData, Inlength, Lfilter, Hfilter,
  152.                                                 filtlen, OutData[2 * i], OutData[(2 * i) + 1]);
  153.         InData = OutData[2 * i];
  154.     }
  155. }
  156.  
  157.  
  158. double DotpEven(double *data, double *filter, char filtlen)
  159. {
  160.     char i;
  161.     double sum;
  162.  
  163.     sum = 0.0;
  164.     for (i = 0; i < filtlen; i += 2)
  165.         sum += *data-- * filter[i]; /* decreasing data pointer is moving */
  166.                                                                 /* backward in time */
  167.     return sum;
  168. }
  169.  
  170.  
  171. double DotpOdd(double *data, double *filter, char filtlen)
  172. {
  173.     char i;
  174.     double sum;
  175.  
  176.     sum = 0.0;
  177.     for (i = 1; i < filtlen; i += 2)
  178.         sum += *data-- * filter[i];    /* decreasing data pointer is moving */
  179.                                                                 /* backward in time */
  180.     return sum;
  181. }
  182.  
  183.  
  184. void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
  185.                     char filtlen, char sum_output, double *output_sequence)
  186.     /* insert zeros between each element of the input sequence and
  187.        convolve with the filter to interpolate the data */
  188. {
  189.     int i;
  190.  
  191.     if (sum_output) /* summation with previous convolution if true */
  192.     {
  193.             /* every other dot product interpolates the data */
  194.         for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
  195.         {
  196.             *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  197.             *output_sequence++ += DotpEven(input_sequence + i + 1, filter, filtlen);
  198.         }
  199.  
  200.         *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  201.     }
  202.     else /* first convolution of pair if false */
  203.     {
  204.             /* every other dot product interpolates the data */
  205.         for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
  206.         {
  207.             *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  208.             *output_sequence++ = DotpEven(input_sequence + i + 1, filter, filtlen);
  209.         }
  210.  
  211.         *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  212.     }
  213. }
  214.  
  215.  
  216. int ReconstructBranches(double *InL, double *InH, int Inlen,
  217.                 double *Lfilter, double *Hfilter, char filtlen, double *Out)
  218.     /* take input data and filters and form two branches of have the
  219.        original length. length of branches is returned */
  220. {
  221.     ConvolveInt2(InL, Inlen, Lfilter, filtlen, 0, Out);
  222.     ConvolveInt2(InH, Inlen, Hfilter, filtlen, 1, Out);
  223.     return Inlen * 2;
  224. }
  225.  
  226.  
  227. void WaveletReconstruction(double **InData, int Inlength, double *Lfilter,
  228.                                 double *Hfilter, char filtlen, char levels, double *OutData)
  229.     /* assumes that the input data has 2 ^ (levels) data points per unit
  230.        interval */
  231. {
  232.     double *Output;
  233.     char i;
  234.  
  235.     Inlength = Inlength >> levels;
  236.     /* destination of all but last branch reconstruction is the next
  237.          higher intermediate approximation */
  238.     for (i = levels - 1; i > 0; i--)
  239.     {
  240.         Output = InData[2 * (i - 1)];
  241.         Inlength = ReconstructBranches(InData[2 * i], InData[(2 * i) + 1],
  242.                                                     Inlength, Lfilter, Hfilter, filtlen, Output);
  243.     }
  244.  
  245.     /* destination of the last branch reconstruction is the output data */
  246.     ReconstructBranches(InData[0], InData[1], Inlength, Lfilter, Hfilter,
  247.                                                                                                             filtlen, OutData);
  248. }
  249.  
  250.  
  251. double CalculateMSE(double *DataSet1, double *DataSet2, int length)
  252. {
  253.     /* calculate mean squared error of output of reconstruction as
  254.          compared to the original input data */
  255.     int i;
  256.     double pointdiff, topsum, botsum;
  257.  
  258.     topsum = botsum = 0.0;
  259.     for (i = 0; i < length; i++)
  260.     {
  261.         pointdiff = DataSet1[i] - DataSet2[i];
  262.         topsum += pointdiff * pointdiff;
  263.         botsum += DataSet1[i] * DataSet1[i];
  264.     }
  265.  
  266.     return topsum / botsum;
  267. }
  268.